Dynamic Regulation Genes at Microtubule Plus Ends: A Novel Class of Glioma Biomarkers

Simple Summary Microtubule plus-end-related genes (MPERGs) encode a group of proteins that specifically aggregate at the microtubule plus ends to play critical biological roles in the cell cycle, cell movement, ciliogenesis, and neuronal development by coordinating microtubule assembly and dynamics; however, the MPERG correlations and their clinical significance in glioma are not fully understood. This study is the first to systematically analyze and define a seven-gene signature (CTTNBP2, KIF18A, NAV1, SLAIN2, SRCIN1, TRIO, and TTBK2) and nomogram model closely associated with clinical factors and the tumor microenvironment as a reliable and independent prognostic biomarker to guide personalized choices of immunotherapy and chemotherapy for glioma patients. Abstract Glioma is the most prevalent and aggressive primary nervous system tumor with an unfavorable prognosis. Microtubule plus-end-related genes (MPERGs) play critical biological roles in the cell cycle, cell movement, ciliogenesis, and neuronal development by coordinating microtubule assembly and dynamics. This research seeks to systematically explore the oncological characteristics of these genes in microtubule-enriched glioma, focusing on developing a novel MPERG-based prognostic signature to improve the prognosis and provide more treatment options for glioma patients. First, we thoroughly analyzed and identified 45 differentially expressed MPERGs in glioma. Based on these genes, glioma patients were well distinguished into two subgroups with survival and tumor microenvironment infiltration differences. Next, we further screened the independent prognostic genes (CTTNBP2, KIF18A, NAV1, SLAIN2, SRCIN1, TRIO, and TTBK2) using 36 prognostic-related differentially expressed MPERGs to construct a signature with risk stratification and prognostic prediction ability. An increased risk score was related to the malignant progression of glioma. Therefore, we also designed a nomogram model containing clinical factors to facilitate the clinical use of the risk signature. The prediction accuracy of the signature and nomogram model was verified using The Cancer Genome Atlas and Chinese Glioma Genome Atlas datasets. Finally, we examined the connection between the signature and tumor microenvironment. The signature positively correlated with tumor microenvironment infiltration, especially immunoinhibitors and the tumor mutation load, and negatively correlated with microsatellite instability and cancer stemness. More importantly, immune checkpoint blockade treatment and drug sensitivity analyses confirmed that this prognostic signature was helpful in anticipating the effect of immunotherapy and chemotherapy. In conclusion, this research is the first study to define and validate an MPERG-based signature closely associated with the tumor microenvironment as a reliable and independent prognostic biomarker to guide personalized choices of immunotherapy and chemotherapy for glioma patients.


Introduction
Glioma derived from neuroglial cells and neurons are the most prevalent as well as aggressive primary intracranial nervous system tumors, accounting for 40-50% of intracranial tumors. The annual incidence of glioma in China is 3-6.4 per 100,000 individuals, accounting for 23.3% of central nervous system (CNS) tumors as well as 78.3% of malignant CNS tumors, and the annual mortality rate is about 30,000 individuals [1,2]. According to the 2016 World Health Organization (WHO) classification, glioma can be divided into four different grades based on the malignancy level. Among them, WHO grade I is benign glioma with a long survival period; WHO grades II and III belong to low-grade glioma (LGG) with a median survival of 5-10 years; and WHO grade IV, also known as glioblastoma (GBM), is a high-grade glioma with high invasive and therapeutic resistance and less than two years' median survival [3]. Genetic and environmental factors are the main risk factors, such as genetic mutation, nitrite-containing food, virus, or bacterial infection, especially high-dose ionizing radiation. Currently, surgical resection plus radiotherapy, chemotherapy, as well as tumor-treating fields are the main treatment methods for glioma. Despite the availability of multiple biomarkers for glioma diagnosis and therapy, the clinical outcome remains poor [4]. Therefore, more sensitive biomarkers and effective therapeutic strategies are critically required to curb glioma progression and improve patients' outcomes.
The tumor microenvironment (TME) primarily contains peripheral blood vessels, immune and stromal cells, as well as non-cellular substances. At different stages of tumor development, different microenvironment components play a leading role, thus affecting the prognosis and treatment of tumors [5]. For instance, stromal components tend to form a tumor-promoting microenvironment, while immune components are initially recruited by tumor cells and form anti-tumor immune microenvironments. However, as tumors grow, the effector immune cells gradually become exhausted and begin to produce a suppressed immune microenvironment to promote tumor growth and invasion. The blood-brain barrier is initially considered a limiting factor for preventing immune cells from entering the CNS. However, with the discovery of lymphatic components in the CNS, the concept of the CNS as an immune-blind area has gradually weakened, rendering immunotherapy a possible and promising treatment option for glioma [6,7].
A microtubule (MT) is one of the most important cytoskeleton components for maintaining cell morphology, and it is also the main molecular target of tumor drug therapy [8][9][10]. As the core regulator of MT dynamics and the interaction between MT and other organelles or structures, MT plus-end binding proteins (MPEBPs) are a group of proteins that specifically aggregate at the plus end of MTs, and their biological functions in various cell types and stages have been extensively studied [11][12][13]. For instance, in mitotic cells, MPEBPs coordinate mitotic spindle dynamics to ensure bipolar orientation and the precise separation of chromosomes, which is critical for maintaining genomic stability. In non-mitotic cells, MPEBPs are mainly responsible for regulating cell polarization, migration, communication, and movement to maintain the normal structure of cells, contributing to the normal occurrence of various physiological processes, including ciliogenesis and neuronal development. Mutation or abnormal regulation of MT plus-end-related genes (MPERGs) leads to a series of MT-related diseases, including ciliopathies and neurodevelopmental and neurodegenerative diseases [12,14]. At present, a growing number of MPERGs have been found to be expressed abnormally in malignant tumors, thus affecting the growth, migration, and invasion of tumors. The dynamic change in the MT plus end also determines the sensitivity of MT-targeted drugs, resulting in different drug responses in patients [15]. Although it is relatively convenient to analyze the function of individual MPERGs in tumors, the systematic study of the multi-gene signature in glioma could better reveal the impact of their overall Biology 2023, 12, 488 3 of 24 relationship on glioma. Additionally, targeted treatment for individual MPERGs is usually insufficient to achieve the expected clinical results. Therefore, the comprehensive study of MPERGs will guide multi-target combined therapy, thus benefiting more glioma patients.
In this study, we used public datasets to systematically analyze MPERG expression, focusing on their clinical significance in glioma. We aimed to develop a novel MPERGbased prognostic signature and evaluate its potential to anticipate the clinical outcome and therapeutic effect in glioma patients. Figure 1 presents a flowchart of the research.
affecting the growth, migration, and invasion of tumors. The dynamic change in the MT plus end also determines the sensitivity of MT-targeted drugs, resulting in different drug responses in patients [15]. Although it is relatively convenient to analyze the function of individual MPERGs in tumors, the systematic study of the multi-gene signature in glioma could better reveal the impact of their overall relationship on glioma. Additionally, targeted treatment for individual MPERGs is usually insufficient to achieve the expected clinical results. Therefore, the comprehensive study of MPERGs will guide multi-target combined therapy, thus benefiting more glioma patients.
In this study, we used public datasets to systematically analyze MPERG expression, focusing on their clinical significance in glioma. We aimed to develop a novel MPERGbased prognostic signature and evaluate its potential to anticipate the clinical outcome and therapeutic effect in glioma patients. Figure 1 presents a flowchart of the research.

Differential Expression Analyses
The distribution of the gene expression was displayed by complex heatmaps using the ComplexHeatmap R package (v 2.2.0) [22]. The comparison of the gene expression in different tissues or groups was displayed using histograms through the Wilcoxon ranksum test. Immunohistochemical (IHC) images of the signature protein expression in healthy and glioma samples were extracted from the Human Protein Atlas (HPA) portal (https://www.protein.org, accessed on 20 June 2022).

Unsupervised Consensus Clustering
Unsupervised consensus clustering of the glioma patients was carried out using the ConsensusClusterPlus R package (v 1.54.0) with 1000 repetitions, and the cumulative distribution function (CDF) curve was employed to gain the appropriate cluster number. A principal component analysis (PCA) was employed to display the distribution differences among the clusters. A KM plot indicating the survival difference was obtained with the survivin R package (v 3.2-10) and displayed with the survminer R package (v 0.4.9) [24][25][26].

Functional Enrichment Analyses
First of all, the MPERGs were directly used for the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional studies with the clusterProfiler R package (v 3.14.3). The cluster-associated genes were identified with the limma R package (v 3.40.2) and visualized using a volcano plot [27]. Then, the cluster-associated genes were further used for a Gene Set Enrichment Analysis (GSEA) containing the h.all.v7.2 (Hallmark) and c2.cp.v7.2 (KEGG) gene sets with the clusterProfiler R package [28,29].

Development and Validation of Risk Signature and Nomogram Model
The independent prognostic genes for constructing the risk signature were screened from 45 DEMPERGs using the univariate, KM, least absolute shrinkage and selection operator (LASSO), and multivariate Cox regression methods in turn with the survival and glmnet R (v 4.1-2) packages [17]. The histograms representing the risk coefficient and concordance index (C-index) were obtained from the Cox regression analyses. The receiver operating characteristic (ROC) curves representing the prognostic value were generated with the pROC (v 1.17.0.1) and timeROC R (v 0.4) packages, while the risk score plot was generated with the ggrisk R package (v 1.3). The correlation among the clusters, risk scores, and survival statuses was displayed by means of a Sankey plot with the ggalluvial R package (v 0.12.3).
The distributions of the risk signature gene expressions and risk scores in diverse clinical statuses were visualized using complex heatmaps. The risk grouping under different clinical factors was compared utilizing a chi-squared test and displayed using a baseline datasheet. To finally screen the independent prognostic factors, uni-and multi-variate Cox regression analyses were carried out and displayed by means of a forest plot. The nomogram model was then developed based on these factors using the rms (v 6.2-0) and survival R packages. Calibration curves and the C-index were used to assess the accuracy of the prediction. The decision curve analysis (DCA) plot for assessing the clinical net benefit was generated with the survival R package and stdca R [30].

TME Infiltration Analyses
The TME infiltration, including the stromal and immune cell infiltration, was calculated via the ssGSEA method in the GSVA package (v 1.34.0), the Estimate package (v 1.0.13), and the xCell method in the immunedeconv package, and visualized using complex heatmaps [31][32][33][34]. The infiltrations in various cells or groups were compared utilizing the Wilcoxon rank-sum test and displayed via histogram. The Spearman correlations between the signature as well as the TME infiltration, tumor mutation burden (TMB), and microsatellite instability (MSI) were examined and visualized using heatmaps [35][36][37]. A stemness analysis was conducted using the OCLR algorithm, and the Spearman correlations between the signature and stemness were visualized using heatmaps [38].

Immunotherapy and Drug Sensitivity Analyses
The Tumor Immune Dysfunction and Exclusion (TIDE, https://tide.dfci.harvard.edu/, accessed on 20 June 2022) portal uses immune escape and immune cell dysfunction to anticipate the immune checkpoint blockade (ICB) treatment response of patients. The TIDE scores in different groups were compared utilizing Welch's t-test and displayed via histograms. High scores indicated stronger immune escape, resulting in poorer ICB efficacy and shorter survival, and vice versa [39,40].
The Gene Set Cancer Analysis (GSCA, http://bioinfo.life.hust.edu.cn/GSCA/#/, accessed on 20 June 2022) portal containing Cancer Therapeutics Response Portal (CTRP) data was employed to anticipate the drug sensitivity according to the connection between the gene expression and half maximal inhibitory concentration (IC50) of small molecules [41]. The correlation heatmap was downloaded from the portal. Positive correlation revealed that increased gene expression was associated with the high IC50 of the drugs, indicating low drug sensitivity and poor therapeutic effect, and vice versa.  [11][12][13]42,43]. Table 1 presents the gene and protein names and MT plus end binding modes of the MPERGs. CAP-Gly domain-containing linker protein 4

Systematic Analyses of MPERG Expression, Genetic Alteration, Correlation, and Interaction in Glioma
Microtubule-associated protein RP/EB family member 1 Autonomous

MAPRE3 (EB3)
Microtubule-associated protein RP/EB family member 3 Transforming acidic coiled-coil-containing protein 3 Unclear To study the correlation and function of these genes in glioma, we systematically analyzed them from the aspects of genetic alteration, DNA methylation, mRNA expression, mRNA correlation, protein-protein interaction, and functional correlation. Firstly, we visualized the genetic alteration of the MPERGs, including the single nucleotide alteration (SNA) and can, from cBioPortal. Only 16% (93/576) of patients had MPERG mutations, and the mutation frequency of individual genes was very low (≤1.7%). The mutation styles mainly included in-frame, missense, splice, and truncating mutations (Figure 2A). To further observe the mutations of other genes, we divided the sample into the MPERG mutant group and the non-mutant group. As shown in the histogram, IDH1, TP53, ATRX, PTEN, EGFR, CIC, TTN, FUBP1, MUC16, and NOTCH were the ten genes with the highest mutation frequency in any group ( Figure 2B; Table S1). Among them, TP53, ATRX, PTEN, EGFR, and NOTCH1 are the 50 most commonly mutated genes in cancer, while IDH1 mutation is a mark event for the early development of glioma, indicating a favorable prognosis [44,45]. However, only the frequency of the CIC mutation had a significant difference between the two groups (p = 2.932 × 10 −4 ). In addition, the overall survival (OS, p = 1.150 × 10 −2 ) and disease-free survival (DFS, 1.484 × 10 −3 ) in the mutant group were shorter than those in the non-mutant group ( Figure 2C). Next, we examined the CNA of the MPERGs. Here, 24% (259/1090) of the patients had MPERG CNAs, and amplification (16%) was more likely to occur than deep deletion (8%) ( Figure 2D). Similarly, the CNA frequency of each gene was very low (<3%) ( Figure 2E; Table S2). Among these patients, CDKN2A, CDKN2B, CDKN2A-DT, MTAP, EGFR, SEC61G, RN7SL151P, DMRTA1, M1R31HG, and IFNE were the ten genes with the highest CNA frequencies, and the frequencies of CDKN2A, CDKN2B, CDKN2A-DT, MTAP, EGFR, SEC61G, and DMRTA1 were lower in the MPERG CNA group than in the non-CNA group ( Figure 2F). Additionally, compared with the non-CNA group, the MPERG CNA group had higher OS (p = 6.097 × 10 −4 ), while no significant variation in the DFS was detected ( Figure 2G). According to our findings, although the genetic alteration of individual MPERGs is very low, the overall alteration of the MPERGs still causes dysfunction, leading to significant prognostic differences in glioma patients.
Previous studies have demonstrated that CNA and DNA methylation are crucial factors in regulating gene expression [46,47]. Our correlation results also showed that most MPERGs positively correlated with the CNAs, while they negatively correlated with DNA methylation ( Figure 3A; Table S3). Therefore, we further evaluated the MPERG expression between the normal and glioma tissues with TCGA and GTEx data from the UCSC Xena database. Based on the Wilcoxon rank-sum test, 45 differentially expressed MPERGs (DEMPERGs) were found, including 36 upregulated and 9 downregulated genes ( Figure 3B,C). Most MPERGs had correlation with each other in glioma ( Figure 3D; Table S4). Additionally, the PPI analysis identified MAPRE1, CKAP5, MAPRE3, CLIP1, KIF11, and KIF2C as hub proteins of interaction (node degree ≥ 10) ( Figure 3E; Table S5).
To further explore the effect of their correlation and interaction on glioma, we conducted a functional enrichment analysis. The GO functional annotations indicated that the MPERGs were primarily located at the MT, cell leading edge, and ruffle, and they were significantly enriched in the cell cycle, epithelial cell migration (ECM), tissue migration, wound healing, and antigen processing presentation of exogenous antigen. Additionally, the KEGG enrichment analysis indicated that they were enhanced in the Hippo signaling pathway, regulation of actin cytoskeleton, Wnt signaling pathway, and signaling pathways regulating the pluripotency of stem cells. It is noteworthy that the last two pathways are also cancer-stemness-related pathways, implying that the MPERGs may participate in cancer-and immune-related pathways, especially cancer-stemness-related pathways ( Figure 3F; Table S6).

Glioma Patients with DEMPERGs Were Well Distinguished into Two Subgroups with Survival and TME Infiltration Differences
Based on the above-mentioned 45 DEMPERGs, the glioma patients were distinguished into two to six subgroups, where k = 3 was the ideal clustering with the most stable k value ( Figures 4A,B and S1A,B). Nevertheless, non-significant variation existed in the OS between clusters 1 and 2 ( Figure S1C). Therefore, we re-clustered the patients according to k = 2 for the subsequent analyses ( Figure 4C). Forty-three of the 45 genes were expressed differently in two clusters ( Figure 4D). The PCA plot confirmed that these two clusters were well distinguished ( Figure 4E), and worse OS in cluster 1 was observed than in cluster 2 ( Figure 4F). Previous studies have demonstrated that CNA and DNA methylation are crucial factors in regulating gene expression [46,47]. Our correlation results also showed that most MPERGs positively correlated with the CNAs, while they negatively correlated with DNA methylation ( Figure 3A; Table S3). Therefore, we further evaluated the MPERG expression between the normal and glioma tissues with TCGA and GTEx data from the UCSC Xena database. Based on the Wilcoxon rank-sum test, 45 differentially expressed MPERGs (DEMPERGs) were found, including 36 upregulated and 9 downregulated genes ( Figure  3B,C). Most MPERGs had correlation with each other in glioma ( Figure 3D; Table S4). Additionally, the PPI analysis identified MAPRE1, CKAP5, MAPRE3, CLIP1, KIF11, and KIF2C as hub proteins of interaction (node degree ≥ 10) ( Figure 3E; Table S5). Subsequently, functional enrichment with GSEA was conducted to evaluate the possible mechanism of survival differences. The result showed that 1040 upregulated and 860 downregulated genes were detected between clusters 1 and 2 (|log 2 fold change| > 1, adjusted p < 0.05) ( Figure 4G). The cluster-associated genes were positively correlated with the cancer-related hallmarks (e.g., epithelial-mesenchymal transition, G2/M checkpoint, apoptosis, angiogenesis, and p53 pathway) and immune-related hallmarks (e.g., IL6-JAK/STAT3 signaling, interferon response, coagulation, and inflammatory response). From the KEGG, they were strongly related to the cancer-related pathways (e.g., cell cycle, ECM receptor interaction, and p53 and JAT-STAT signaling pathways) and immune-related processes (e.g., antigen processing and presentation, and cytokine-cytokine receptor interaction), and negatively related to the MAPK, Wnt, and ERBB signaling pathways. Among them, the JAK-STAT, MAPK, and Wnt signaling pathways are also cancer-stemness-related pathways ( Figure 4H; Table S7). These findings indicated that besides the cell cycle and cell migration, the MPERG-associated genes also participate in the cancer-and immuneassociated pathways.
Considering the close relationship between the MPERG-associated genes and abovementioned pathways, we investigated the TME infiltration between the two clusters. Based on the ssGSEA method, the T, cytotoxic T lymphocyte (CTL), CD4 T helper (Th) 1/2/17, natural killer (NK), NK CD56 dim cells, immature dendritic cells (iDCs), activated DCs (aDCs), macrophages, neutrophils, and eosinophils were strongly infiltrated in the cluster 1 patients. In contrast, the cluster 2 patients were mainly infiltrated in the CD8 T, CD4 follicle helper T (Tfh), CD4 regulatory T (Treg), gamma delta T (Tgd), central memory T (Tcm), effector memory T (Tem), NK CD56 bright cells, and plasmacytoid DCs (pDCs) ( Figure 5A,B). In general, the cluster 1 patients had markedly higher TME scores than the cluster 2 patients, including the stromal, immune, and ESTIMATE scores, as evaluated using the ESTIMATE method ( Figure 5A,C). Next, we used the xCell algorithm to verify this result. The patients in cluster 1 showed high infiltrations in the CD8 T naive/Tcm/Tem, CD4 Tm/Tem/Th1/Th2 cells, myeloid DC (mDC), pDC, monocyte, macrophage, macrophage M1, eosinophil, common lymphoid progenitor, and endothelial cells. The cluster 2 patients had high infiltrations in the B cell plasma, B cell (class-switched memory), CD4 Treg, NK, mast cells, T cell NK, activated myeloid DCs (amDCs), neutrophils, and granulocytemonocyte progenitor cells ( Figure 5D,E). The TME scores were the same as the ESTIMATE algorithm scores (Figure 5D,F). Next, we investigated the immune-related gene expression in the two clusters. The expression of the major histocompatibility complexes and most immunostimulators in cluster 1 was higher than that in cluster 2. Surprisingly, most of the immunoinhibitor expression in cluster 1 was also higher, implying that there is a relatively complex immune microenvironment in MPERG-related glioma, and their relationship with prognosis and treatment needs further analysis ( Figure 5G-I). In conclusion, the TME infiltration was significantly different between the two subgroups of glioma patients.

A Seven-Gene Prognostic Signature Was Constructed and Validated in Glioma
To screen the predictive factors from the 45 DEMPERGs, we first conducted univariate and KM Cox regression analyses. A total of 11 risk genes and 25 protective genes were identified as having predictive value for the OS, disease-specific survival, as well as progression-free interval (Figures 6A and S2). Therefore, a LASSO Cox regression analysis was performed on these 36 genes. A total of 15 genes were selected after taking the variables with minimized lambda (Figure 6B,C). The multivariate Cox regression analysis of these 15 genes indicated that 7 of them were independent prognostic genes ( Figure 6D). Thus, we established a prognostic signature based on these seven gene. Figure 6E). TTBK2, NAV1, and CTTNBP2 were considered as protective factors, while SRCIN1, TRIO, KIF18A, and SLAIN2 were risk factors for glioma. Compared with a single gene, the sevengene signature had the highest C-index in the Cox regression analyses (0.828) and area under the curve (AUC) (0.814) in the ROC curve, indicating that it is the most appropriate factor for predicting glioma patients' clinical outcome ( Figure 6F,G). The Sankey plot showed that the patients in cluster 1 had high risk scores and dead status, while the patients in cluster 2 had the opposite ( Figure 6H).
We further used the TCGA and CGGA datasets as training and validation sets to assess and validate the risk signature. All the signature genes were expressed differentially between the two groups ( Figure 6I,M). Moreover, with the risk score elevated, the survival time became shorter as well as closer to death ( Figure 6J,N). The time-dependent ROC analysis revealed that the 2-, 4-, and 6-year OS AUCs in the TCGA were 0.905, 0.881, and 0.871, while the AUCs in the CGGA were 0.832, 0.853, and 0.857, respectively ( Figure 6K,O). The KM plot confirmed that the high-risk patients showed a worse OS in comparison with the low-risk patients ( Figure 6L,P). In addition, we further verified the signature protein expression in the normal and glioma samples using the HPA database. Compared with the normal samples, except for SLAIN2, the CTTNBP2, KIF18A, NAV1, and TRIO protein expression in glioma was elevated, while SRCIN1 and TTBK2 were decreased ( Figure 6Q). Referring to Figure 3B, the signature gene was basically consistent with their protein expression. These findings confirmed that the risk signature we constructed has very high prediction accuracy.

A Predictive Nomogram Model Was Established and Verified to Facilitate the Clinical Application of the Risk Signature in Glioma
Due to the strong correlation between the risk score and prognosis, we also incorporated clinical factors into the prognosis model. According to the TCGA dataset, an increased risk score was highly related to the malignant progression of glioma, including IDH-WT status, 1p/19q non-codeletion, older age, higher WHO grade, more inclined to form glioblastoma, and worse response to treatment (Table 2; Figure 7A). A similar result was also obtained from the CGGA datasets (Table 3; Figure 7G). To determine whether the signature and clinical features are independent prognostic factors in glioma, we conducted Cox regression analyses and observed that the risk score, PTO, as well as age had independent prognostic value in the TCGA dataset ( Figure 7B), while the CGGA datasets presented more factors, including the WHO grade, 1p/19q codeletion, PRS type, age, and risk score ( Figure 7H). Therefore, based on these independent prognostic factors, we finally established a nomograph model for predicting the 2-, 4-, and 6-year OS ( Figure 7C,I). The results showed that the calibration curve matched the ideal diagonal line well ( Figure 7D,J), and the C-index of the nomogram was also better than a single clinical factor or risk signature (0.884 for TCGA; 0.792 for CGGA) ( Figure 7E,K). In addition, the DCA showed that the nomogram was the best prognostic indicator ( Figure 7F,L). These results demonstrated that our constructed and validated nomogram including clinical factors is the best predictor of glioma prognosis.

Signature-Associated Glioma Patients Were Predicted to Benefit from Immunotherapy and Chemotherapy
Moreover, we continued our investigation into the connection between the risk signature and TME. The ssGSEA method revealed that the infiltrations of the T, CTL, CD4 Th1/2/17, NK, and NK CD56 dim cells, iDCs, aDCs, macrophages, neutrophils, and eosinophils were strongly connected with the risk scores, whereas the CD8 T, CD4 Tfh, CD4 Treg, Tgd, Tcm, Tem, and NK CD56 bright cells, DC and pDC were inversely correlated with them. Overall, the TME scores evaluated using the ESTIMATE algorithm were positively related to the risk scores ( Figures 8A,B and S3A,B; Table S8). In addition, similar results were obtained using the xCell algorithm ( Figures 8C,D and S3C,D; Table S9).
Subsequently, we continued to explore the response of the glioma patients to immunotherapy and chemotherapy. Immunotherapy is known to be influenced by immunoinhibitor expressions, TMB and MSI scores and other factors [36,37]. Among the 24 immunoinhibitors, 21 of them were strongly related to the risk score ( Figures 8E and S3E; Table S10). High TMB and MSI scores indicate that tumor cells have higher antigenicity and more neoantigens, thus having stronger immune infiltration. The outcomes indicated a positive correlation between the risk score and TMB as well as a negative relationship to the MSI (Figures 8F and S3F; Table S11). These results strongly suggested that immunotherapy is effective for glioma. Therefore, TIDE was utilized to anticipate the efficacy of ICB therapy. The result showed that the TIDE score was lesser in the low-compared to the high-risk patients, indicating a better immunotherapy response and treatment outcome for the low-risk patients ( Figure 8G,H).
Cancer stem cells are mainly accountable for the maintenance and growth of cancer cells, which are highly correlated with tumor metastasis, recurrence, as well as therapy resistance [48]. The functional enrichment analysis suggested that the MPERGs and MEPRG-associated genes may be involved in cancer stemness. Therefore, we introduced the OCLR algorithm to assess the variations of cancer stemness between both groups. The results presented a negative relationship between the risk and stem score, indicating that the low-risk patients had stronger stem cell characteristics and differentiation potential than the high-risk patients ( Figure 8I).
Gene expression may affect the clinical response of patients to drug therapy; hence, these genes represent possible drug screening targets. We studied the connection between the signature expression and drug sensitivity utilizing the CTRP dataset from the GSCA site. The findings revealed that TTBK2 and SRCIN were inversely associated with the IC50 of the majority of the drugs, while TRIO, NAV1, and CTTNBP2 positively correlated with them, implying that TTBK2 and SRCIN increased expression or TRIO, NAV1, and CTTNBP2 reduced expression could enhance the sensitivity of patients to these drugs ( Figure 8J). Therefore, these genes could be new chemotherapeutic targets.

Discussion
Despite treatment advances, glioma remains a fatal disease. Therefore, new therapeutic targets and methods are urgently needed. In recent years, more and more MPERGs have been found to be aberrantly expressed in several cancers, which could be connected to tumor occurrence, development, and prognosis. In our research, we first thoroughly analyzed 46 MPERGs regulated by CNA and DNA methylation and found that 45 MPERGs were differentially expressed in MT highly enriched glioma. Glioma patients with 45 DEM-PERGs could be distinguishably divided into two subgroups with significant prognostic and functional differences. Next, we used univariate, KM, LASSO, and multivariate Cox regression analyses in turn to finally screen seven independent prognostic genes and then constructed the prognostic signature based on the MPERGs for the first time.
The biological functions of these seven signature genes have been extensively studied, although their correlation with cancer, especially glioma, is still limited. For example, CTTNBP2 is an autism-related gene mainly expressed in neurons and highly enriched in the dendritic spine. It is involved in stabilizing microtubules and controlling dendritic spine formation [49]. Upregulation of CTTNBP2 promotes the cancerous growth of the ovary [50]. KIF18A is an MT depolymerizing kinesin that ensures chromosome stability during mitosis by inhibiting MT dynamics without disrupting their stability [51]. KIF18A upregulation contributes to the malignant progression and unfavorable prognosis of various types of cancer, including breast, hepatocellular, clear cell renal, prostate, esophageal carcinomas, and lung adenocarcinoma, especially glioma [52][53][54][55][56][57][58][59][60][61]. SLAIN2 promotes the nucleation and elongation of cytoplasmic MT and suppresses their catastrophes, which is necessary for maintaining the normal structure of the cytoskeleton [62]. Increased SLAIN2 expression encourages the invasion of mesenchymal cells in three-dimensional culture and mouse cancer models. In addition, it promotes MT migration and increases the invasion as well as metastasis of colon cancer, leading to a poor prognosis [63,64]. SRCIN1 plays a role in synaptic maintenance and neurotransmitter release. Moreover, as an anti-tumor molecule, its downregulation stimulates the spread, migration, and invasion of lung, liver, breast, and gastric cancers and neuroblastoma, thereby affecting their prognosis and treatment [65][66][67][68][69]. TTBK2 inhibits the activity of the microtubule depolymerase of KIF2A through phosphorylation, thus promoting the migration of cervical cancer cells. Additionally, it also has a crucial function in cilia formation [70,71]. TTBK2 downregulation sensitizes melanoma cells and renal cells to sunitinib [72]. NAV1 is essential for neurogenesis and participates in neuronal migration [73]. TRIO promotes actin remodeling, cell growth and migration, which is necessary for maintaining proper axonal growth and MT network stability [74][75][76]. However, the correlation between NAV1, TRIO, and cancer has not been well investigated. Our study found that compared with normal tissues, CTTNBP2, KIF18A, NAV1, SLAIN2, and TRIO were upregulated in glioma, while SRCIN1 and TTBK2 were downregulated. Moreover, TTBK2, NAV1, and CTTNBP2 were protective factors for glioma, while SRCIN1, TRIO, KIF18A, and SLAIN2 were risk factors. In general, the signature composed of these seven genes is a risk factor for glioma development, and with an increase in the risk score, glioma was more inclined to malignant progression. Therefore, we also designed a nomogram model containing clinical features and risk scores. Utilizing the TCGA and CGGA datasets, we further confirmed that our construction is the best prognostic predictor for glioma patients.
In addition to regulating MT plus end dynamics to affect glioma, this research also studied the connection between the clusters, risk groups, and TME infiltration. Consistent with the cluster results, the risk scores were strongly connected to the TME scores, especially the immune score, indicating that they may form an anti-tumor immune microenvironment. Nonetheless, a growing number of investigations have pointed out that different immune cells exhibit different immune functions, which may promote or inhibit the antitumor immune response. CD4 Th2, CD4 Treg, and macrophage M2 are known to be representative immunosuppressive cells that promote cancer growth [77][78][79]; hence, we further studied these cells. The risk score was strongly connected to the CD4 Th2 cell infiltration and negatively connected to the CD4 Treg cell infiltration by the ssGSEA algorithm. In addition to these two consistent correlations, the xCell algorithm also revealed that the risk score was strongly connected to macrophage M2, implying that the anti-tumor immune microenvironment is accompanied by the tumor-promoting immune microenvironment. In general, we speculated that in the signature-associated glioma patients, with the malignant progression of the disease, the increased risk score leading to an unfavorable prognosis is related to the suppressive immune microenvironment.
As a novel and promising treatment method, immunotherapy for glioma has attracted increasing attention. Our results indicated that the risk signature was strongly connected to the immunoinhibitors and TMB, and negatively connected to the MSI in glioma. In addition, the ICB therapy prediction confirmed that the treatment effect was inferior in the high-risk patients compared to those at low risk. Combined with the clinical factors, these findings revealed that ICB therapy may be more efficient for LGG patients than GBM patients. Therefore, according to the current results, our signature could be utilized as a reliable biomarker for immunotherapy in glioma patients at different stages of progression. Additionally, we found several associations between the IC50 of chemical drugs and signature gene expression, representing another potential molecular target for glioma chemotherapy that may help to predict chemotherapy resistance.
Inevitably, this research also has certain drawbacks. Although the information was gathered from several public databases, additional experimental verification is needed. In addition, the data showed some unexpected results, such as the negative correlation between the cancer stemness and risk score. However, the progression, prognosis and treatment of glioma are inherently related to multiple factors, and the contribution of a single factor to the above-mentioned process is often limited. Therefore, further research in larger-scale and multicenter clinical cohorts is needed to overcome these inconsistent results.

Conclusions
In summary, this work represents the first to construct and validate an MPERG-based prognostic signature associated with the TME as a reliable and independent prognostic biomarker to provide personalized choices of immunotherapy and chemotherapy for glioma patients.